
#ifndef InputFile_H
#define InputFile_H

#include <vector>
#include <fstream>
#include <string>
#include <complex>
#include <cctype>

#include <fitsio.h>

#include "dbg.h"
#include "ConfigFile.h"
#include "Angle.h"

// Helper function to convert ASCII line into vector of values
static bool GetValues(std::istream& is, std::vector<double>& values,
                      char delim, const std::string comment_markers)
{
    // Skip if this is a comment.
    char first = is.peek();
    // Skip initial whitespace
    while (std::isspace(first)) { is.get(first); first = is.peek(); }
    bool skip = false;
    if (comment_markers.find(first) != std::string::npos) skip = true;
    if (skip) {
        is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
        return GetValues(is,values,delim,comment_markers);
    }

    const int nval = values.size();
    if (delim == '\0') {
        double temp;
        for(int i=0;i<nval;++i) {
            is >> temp;
            values[i] = temp;
        }
    } else {
        ConvertibleString temp;
        for(int i=0;i<nval;++i) {
            getline(is,temp,delim);
            values[i] = temp;
        }
    }
    if (!is) return false;
    is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
    return true;
}

static angle::AngleUnit ConvertUnitStr(const std::string& str)
{
    if (str.find("rad") == 0) return angle::radians;
    else if (str.find("deg") == 0) return angle::degrees;
    else if (str.find("hour") == 0) return angle::hours;
    else if (str.find("arcmin") == 0) return angle::arcmin;
    else if (str.find("arcsec") == 0) return angle::arcsec;
    else myerror("Invalid angle unit: ",str.c_str());
    // Won't get here, but this avoids a warning.
    return angle::radians;
}

static void fitserror(int status)
{
    char msg[80];
    std::cerr<<"FITS error found: \n";
    fits_get_errstatus(status,msg);
    std::cerr<<"  "<<msg<<std::endl;
    while (fits_read_errmsg(msg))
        std::cerr<<msg<<std::endl;
    myerror("Fits error");
}

class InputFile
{
public:
    InputFile(const std::string& file_name, const ConfigFile& params) : _use_ra_dec(false)
    {
        if (params.keyExists("file_type")) {
            // If file_type is specified, use that.
            std::string file_type = params["file_type"];
            if (file_type == "ASCII" || file_type == "ascii")
                readAscii(file_name,params);
            else if (file_type == "FITS" || file_type == "fits")
                readFits(file_name,params);
            else 
                myerror("Unknown file_type ",file_type.c_str());
        } else {
            // Otherwise, try to get the type from the file extension.
            if (file_name.find(".fit") != std::string::npos)
                readFits(file_name,params);
            else
                readAscii(file_name,params);
        }
        dbg<<"Input file "<<file_name<<": ngal = "<<_pos.size()<<std::endl;
        xdbg<<"shear.size = "<<_shear.size()<<std::endl;
        xdbg<<"weight.size = "<<_weight.size()<<std::endl;

        // Check if we need to flip g1 or g2
        bool flip1 = params.read<bool>("flip_g1",false);
        bool flip2 = params.read<bool>("flip_g2",false);
        if (flip1 && flip2) {
            dbg<<"flipping sign of both g1 and g2\n";
            for(size_t i=0; i<_shear.size(); ++i) _shear[i] = -_shear[i];
        } else if (flip2) {
            dbg<<"flipping sign of g2\n";
            for(size_t i=0; i<_shear.size(); ++i) _shear[i] = conj(_shear[i]);
        } else if (flip1) {
            dbg<<"flipping sign of g1\n";
            for(size_t i=0; i<_shear.size(); ++i) _shear[i] = -conj(_shear[i]);
        }

        if (_use_ra_dec && params.read("project",false)) {
            // Project the position onto a tangent plane
            double ra_cen, dec_cen;
            if (params.keyExists("project_ra")) {
                ra_cen = params["project_ra"];
                std::string ra_unit_str = params.read("ra_units",std::string("arcsec"));
                ra_cen = (ra_cen * ConvertUnitStr(ra_unit_str)) / angle::radians;
            } else {
                ra_cen = 0.;
                for (size_t i=0;i<_pos.size();++i) ra_cen += _pos[i].getX();
                ra_cen /= _pos.size();
            }
            if (params.keyExists("project_dec")) {
                dec_cen = params["project_dec"];
                std::string dec_unit_str = params.read("dec_units",std::string("arcsec"));
                dec_cen = (dec_cen * ConvertUnitStr(dec_unit_str)) / angle::radians;
            } else {
                dec_cen = 0.;
                for (size_t i=0;i<_pos.size();++i) dec_cen += _pos[i].getY();
                dec_cen /= _pos.size();
            }
            project(ra_cen,dec_cen);
        }
    }

    void readAscii(const std::string& file_name, const ConfigFile& params)
    {
        std::ifstream fin(file_name.c_str());

        int max_col = 0;

        int x_col = params.read<int>("x_col",0);
        int y_col = params.read<int>("y_col",0);
        int ra_col = params.read<int>("ra_col",0);
        int dec_col = params.read<int>("dec_col",0);
        double x_units=1., y_units=1., ra_units=1., dec_units=1.;

        if (x_col > 0 || y_col > 0) {
            std::string x_unit_str = params.read("x_units",std::string("arcsec"));
            x_units = (1. * ConvertUnitStr(x_unit_str)) / angle::radians;
            std::string y_unit_str = params.read("y_units",std::string("arcsec"));
            y_units = (1. * ConvertUnitStr(y_unit_str)) / angle::radians;
            if (x_col <= 0) myerror("x_col missing or invalid");
            if (y_col <= 0) myerror("y_col missing or invalid");
            if (ra_col > 0 || dec_col > 0) 
                myerror("ra/dec cols are not allowed in conjunction with x/y cols");
            if (x_col > max_col) max_col = x_col;
            if (y_col > max_col) max_col = y_col;
        } else if (ra_col > 0 || dec_col > 0) {
            std::string ra_unit_str = params.read("ra_units",std::string("arcsec"));
            ra_units = (1. * ConvertUnitStr(ra_unit_str)) / angle::radians;
            std::string dec_unit_str = params.read("dec_units",std::string("arcsec"));
            dec_units = (1. * ConvertUnitStr(dec_unit_str)) / angle::radians;
            if (ra_col <= 0) myerror("ra_col missing or invalid");
            if (dec_col <= 0) myerror("dec_col missing or invalid");
            _use_ra_dec = true;
            if (ra_col > max_col) max_col = ra_col;
            if (dec_col > max_col) max_col = dec_col;
        } else {
            myerror("No valid position columns specified");
        }

        int e1_col = params.read<int>("e1_col",0);
        int e2_col = params.read<int>("e2_col",0);
        int g1_col = params.read<int>("g1_col",0);
        int g2_col = params.read<int>("g2_col",0);

        if (g1_col > 0 || g2_col > 0) {
            if (g1_col <= 0) myerror("g1_col missing or invalid");
            if (g2_col <= 0) myerror("g2_col missing or invalid");
            if (e1_col > 0 || e2_col > 0) 
                myerror("e1/e2 cols are not allowed in conjunction with g1/g2 cols");
            if (g1_col > max_col) max_col = g1_col;
            if (g2_col > max_col) max_col = g2_col;
        } else if (e1_col > 0 || e2_col > 0) {
            if (e1_col <= 0) myerror("e1_col missing or invalid");
            if (e2_col <= 0) myerror("e2_col missing or invalid");
            if (e1_col > max_col) max_col = e1_col;
            if (e2_col > max_col) max_col = e2_col;
        }

        int w_col = params.read<int>("w_col",0);
        if (w_col > max_col) max_col = w_col;

        // Set up delimiter and allowed comment markers
        char delim = params.read("delimiter",'\0');
        std::vector<char> cm_vec = params.read("comment_marker",std::vector<char>(1,'#'));
        std::string comment_markers(1,cm_vec[0]);
        for (size_t i=1;i<cm_vec.size();++i) comment_markers.push_back(cm_vec[i]);

        std::vector<double> values(max_col);
        while (GetValues(fin,values,delim,comment_markers)) {

            // Position
            if (x_col > 0) {
                Assert(x_col <= int(values.size()));
                Assert(y_col <= int(values.size()));
                double x = values[x_col-1] * x_units;
                double y = values[y_col-1] * y_units;
                _pos.push_back(Position<Flat>(x,y));
            } else {
                Assert(ra_col <= int(values.size()));
                Assert(dec_col <= int(values.size()));
                double ra = values[ra_col-1] * ra_units;
                double dec = values[dec_col-1] * dec_units;
                _pos.push_back(Position<Flat>(ra,dec));
            }

            // Shear
            if (g1_col > 0) {
                Assert(g1_col <= int(values.size()));
                Assert(g2_col <= int(values.size()));
                double g1 = values[g1_col-1];
                double g2 = values[g2_col-1];
                _shear.push_back(std::complex<double>(g1,g2));
            } else {
                Assert(e1_col <= int(values.size()));
                Assert(e2_col <= int(values.size()));
                double e1 = values[e1_col-1];
                double e2 = values[e2_col-1];
                double esq = e1*e1+e2*e2;
                // e = (1-q^2)/(1+q^2)
                // g = (1-q)/(1+q)
                // e+eq^2 = 1-q^2
                // q^2 = (1-e)/(1+e)
                // scale = g/e
                //       = (after some manipulations...)
                //       = 1/(1+sqrt(1-e^2))
                double scale = 1./(1.+sqrt(1.-esq));
                _shear.push_back(std::complex<double>(e1*scale,e2*scale));
            }

            // Weight
            if (w_col > 0) {
                Assert(w_col <= int(values.size()));
                double w = values[w_col-1];
                _weight.push_back(w);
            }
        }
    }

    void readFits(const std::string& file_name, const ConfigFile& params)
    {
        fitsfile* fptr;
        int status=0;
        fits_open_table(&fptr, file_name.c_str(), READONLY, &status);
        if (status) fitserror(status);

        long nrows;
        fits_get_num_rows(fptr,&nrows,&status);
        if (status) fitserror(status);

        std::string x_col = params.read<std::string>("x_col","");
        std::string y_col = params.read<std::string>("y_col","");
        std::string ra_col = params.read<std::string>("ra_col","");
        std::string dec_col = params.read<std::string>("dec_col","");
        double x_units=1., y_units=1., ra_units=1., dec_units=1.;

        if (x_col != "" || y_col != "") {
            std::string x_unit_str = params.read("x_units",std::string("arcsec"));
            x_units = (1. * ConvertUnitStr(x_unit_str)) / angle::radians;
            std::string y_unit_str = params.read("y_units",std::string("arcsec"));
            y_units = (1. * ConvertUnitStr(y_unit_str)) / angle::radians;
            if (x_col == "") myerror("x_col missing or invalid");
            if (y_col == "") myerror("y_col missing or invalid");
            if (ra_col != "" || dec_col != "") 
                myerror("ra/dec cols are not allowed in conjunction with x/y cols");

            int x_colnum, y_colnum, anynul;
            fits_get_colnum(fptr,CASEINSEN,(char*)x_col.c_str(),&x_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> x_vec(nrows);
            fits_read_col(fptr,TDOUBLE,x_colnum,1,1,nrows,0,&x_vec[0],&anynul,&status);

            fits_get_colnum(fptr,CASEINSEN,(char*)y_col.c_str(),&y_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> y_vec(nrows);
            fits_read_col(fptr,TDOUBLE,y_colnum,1,1,nrows,0,&y_vec[0],&anynul,&status);

            _pos.resize(nrows);
            for(int i=0;i<nrows;++i) {
                _pos[i] = Position<Flat>(x_vec[i]*x_units,y_vec[i]*y_units);
            }

        } else if (ra_col != "" || dec_col != "") {
            std::string ra_unit_str = params.read("ra_units",std::string("arcsec"));
            ra_units = (1. * ConvertUnitStr(ra_unit_str)) / angle::radians;
            std::string dec_unit_str = params.read("dec_units",std::string("arcsec"));
            dec_units = (1. * ConvertUnitStr(dec_unit_str)) / angle::radians;
            if (ra_col == "") myerror("ra_col missing or invalid");
            if (dec_col == "") myerror("dec_col missing or invalid");
            _use_ra_dec = true;

            int ra_colnum, dec_colnum, anynul;
            fits_get_colnum(fptr,CASEINSEN,(char*)ra_col.c_str(),&ra_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> ra_vec(nrows);
            fits_read_col(fptr,TDOUBLE,ra_colnum,1,1,nrows,0,&ra_vec[0],&anynul,&status);

            fits_get_colnum(fptr,CASEINSEN,(char*)dec_col.c_str(),&dec_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> dec_vec(nrows);
            fits_read_col(fptr,TDOUBLE,dec_colnum,1,1,nrows,0,&dec_vec[0],&anynul,&status);

            _pos.resize(nrows);
            for(int i=0;i<nrows;++i) {
                _pos[i] = Position<Flat>(ra_vec[i]*ra_units,dec_vec[i]*dec_units);
            }

        } else {
            myerror("No valid position columns specified");
        }

        std::string e1_col = params.read<std::string>("e1_col","");
        std::string e2_col = params.read<std::string>("e2_col","");
        std::string g1_col = params.read<std::string>("g1_col","");
        std::string g2_col = params.read<std::string>("g2_col","");

        if (g1_col != "" || g2_col != "") {
            if (g1_col == "") myerror("g1_col missing or invalid");
            if (g2_col == "") myerror("g2_col missing or invalid");
            if (e1_col != "" || e2_col != "") 
                myerror("e1/e2 cols are not allowed in conjunction with g1/g2 cols");

            int g1_colnum, g2_colnum, anynul;
            fits_get_colnum(fptr,CASEINSEN,(char*)g1_col.c_str(),&g1_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> g1_vec(nrows);
            fits_read_col(fptr,TDOUBLE,g1_colnum,1,1,nrows,0,&g1_vec[0],&anynul,&status);

            fits_get_colnum(fptr,CASEINSEN,(char*)g2_col.c_str(),&g2_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> g2_vec(nrows);
            fits_read_col(fptr,TDOUBLE,g2_colnum,1,1,nrows,0,&g2_vec[0],&anynul,&status);

            _shear.resize(nrows);
            for(int i=0;i<nrows;++i) {
                _shear[i] = std::complex<double>(g1_vec[i],g2_vec[i]);
            }

        } else if (e1_col != "" || e2_col != "") {
            if (e1_col == "") myerror("e1_col missing or invalid");
            if (e2_col == "") myerror("e2_col missing or invalid");

            int e1_colnum, e2_colnum, anynul;
            fits_get_colnum(fptr,CASEINSEN,(char*)e1_col.c_str(),&e1_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> e1_vec(nrows);
            fits_read_col(fptr,TDOUBLE,e1_colnum,1,1,nrows,0,&e1_vec[0],&anynul,&status);

            fits_get_colnum(fptr,CASEINSEN,(char*)e2_col.c_str(),&e2_colnum,&status);
            if (status) fitserror(status);
            std::vector<double> e2_vec(nrows);
            fits_read_col(fptr,TDOUBLE,e2_colnum,1,1,nrows,0,&e2_vec[0],&anynul,&status);

            _shear.resize(nrows);
            for(int i=0;i<nrows;++i) {
                double esq = e1_vec[i]*e1_vec[i]+e2_vec[i]*e2_vec[i];
                double scale = 1./(1.+sqrt(1.-esq));
                _shear[i] = std::complex<double>(e1_vec[i]*scale,e2_vec[i]*scale);
            }
        }

        std::string w_col = params.read<std::string>("w_col","");
        if (w_col != "") {
            int w_colnum, anynul;
            fits_get_colnum(fptr,CASEINSEN,(char*)w_col.c_str(),&w_colnum,&status);
            if (status) fitserror(status);
            _weight.resize(nrows);
            fits_read_col(fptr,TDOUBLE,w_colnum,1,1,nrows,0,&_weight[0],&anynul,&status);
        }
    }

    void project(double ra0, double dec0)
    {
        dbg<<"Projecting points relative to (RA,Dec) = ("<<ra0<<','<<dec0<<") rad";
        dbg<<" = ("<<ra0*180./M_PI<<','<<dec0*180./M_PI<<") dec\n";

        // I use a stereographic projection, which preserves angles, but not distances
        // (you can't do both).  The distance error increases with distance from the 
        // projection point of course.
        // The equations are given at:
        //     http://mathworld.wolfram.com/StereographicProjection.html
        // x = k cos(phi) sin(lam-lam0)
        // y = k ( cos(phi0) sin(phi) - sin(phi0) cos(phi) cos(lam-lam0) )
        // k = 2 ( 1 + sin(phi0) sin(phi) - cos(phi0) cos(phi) cos(lam-lam0) )^-1
        //
        // In our case, lam = ra, phi = dec
        double cosdec0 = cos(dec0);
        double sindec0 = sin(dec0);
        for (size_t i=0; i<_pos.size(); ++i) {
            double ra = _pos[i].getX();
            double dec = _pos[i].getY();
            xdbg<<"Project ("<<ra*180./M_PI<<','<<dec*180./M_PI<<") relative to ("<<
                ra0*180./M_PI<<','<<dec0*180./M_PI<<")\n";

            double cosdec = cos(dec);
            double sindec = sin(dec);
            double cosdra = cos(ra - ra0);
            // Note: - sign here is to make +x correspond to -ra,
            //       so x increases for decreasing ra.
            //       East is left on the sky!
            double sindra = -sin(ra - ra0);
            xxdbg<<"cosdec, sindec = "<<cosdec<<','<<sindec<<std::endl;
            xxdbg<<"cosdra, sindra = "<<cosdra<<','<<sindra<<std::endl;

            double k = 2. / (1. + sindec0 * sindec + cosdec0 * cosdec * cosdra );
            xxdbg<<"k = "<<k<<std::endl;
            double x = k * cosdec * sindra;
            double y = k * ( cosdec0 * sindec - sindec0 * cosdec * cosdra );
            xdbg<<"x,y -> "<<x<<','<<y<<std::endl;

            _pos[i] = Position<Flat>(x,y);

            // Now project the ellipticities if necessary:
            if (_shear.size() > 0) {
                Assert(i < _shear.size());
                // For the projection, we use a spherical triangle with A = the point, 
                // B = the projection center, and C = the north pole.
                // The known values are:
                //   a = d(B,C) = Pi/2 - dec0
                //   b = d(A,C) = Pi/2 - dec
                //   C = ra - ra0
                // The value by which we need to rotate the shear is Pi - (A + B)
                // cos(Pi-(A+B)) = -cos(A+B) = -cosA cosB + sinA sinB
                // sin(Pi-(A+B)) = sin(A+B) = sinA cosB + cosA sinB
                //
                // cosc = cosa cosb + sina sinb cosC
                // cosc = sindec0 * sindec + cosdec0 * cosdec * cosdra;
                // cosa = cosb cosc + sinb sinc cosA
                // cosA = (sindec0 - sindec cosc) / (cosdec sinc)
                //      = (sindec0 - sindec (sindec0 sindec + cosdec0 cosdec cosdra) ) /
                //              (cosdec sinc)
                //      = (sindec0 cosdec - sindec cosdec0 cosdra ) / sinc
                // sinA / sina = sinC / sinc
                // sinA = cosdec0 * sindra / sinc
                //
                // Stable solution in case sinc (or cosdec) is small is to calculate denominators,
                // and then normalize using sinA^2 + cosA^2 = 1
                double cosA = sindec0 * cosdec - sindec * cosdec0 * cosdra;
                double sinA = cosdec0 * sindra;
                double normA = sqrt(cosA*cosA + sinA*sinA);
                cosA /= normA;
                sinA /= normA;
                xxdbg<<"A = "<<atan2(sinA,cosA)<<"  = "<<atan2(sinA,cosA)*180/M_PI<<" degrees\n";

                // cosb = cosa cosc + sina sinc cosB
                // cosB = (sindec - sindec0 cosc) / (cosdec0 sinc)
                //      = (sindec - sindec0 (sindec0 sindec + cosdec0 cosdec cosdra) ) / 
                //              (cosdec0 sinc)
                //      = (sindec cosdec0 - sindec0 cosdec cosdra ) / sinc
                // sinB / sinb = sinC / sinc
                // sinB = cosdec * sindra / sinc
                double cosB = sindec * cosdec0 - sindec0 * cosdec * cosdra;
                double sinB = cosdec * sindra;
                double normB = sqrt(cosB*cosB + sinB*sinB);
                cosB /= normB;
                sinB /= normB;
                xxdbg<<"B = "<<atan2(sinB,cosB)<<"  = "<<atan2(sinB,cosB)*180/M_PI<<" degrees\n";

                double cosbeta = -cosA * cosB + sinA * sinB;
                double sinbeta = sinA * cosB + cosA * sinB;
                xxdbg<<"beta = "<<atan2(sinbeta,cosbeta)*180/M_PI<<std::endl;
                std::complex<double> expibeta(cosbeta,sinbeta);
                xxdbg<<"expibeta = "<<expibeta<<std::endl;
                std::complex<double> exp2ibeta = expibeta * expibeta;
                xdbg<<"shear = "<<_shear[i]<<" * "<<exp2ibeta<<" = "<<_shear[i]*exp2ibeta<<std::endl;
                _shear[i] *= exp2ibeta;
            }
        }
        _use_ra_dec = false;
    }

    template <int M>
    void buildCellData(double minsep, double binsize,
                       std::vector<boost::shared_ptr<CellData<NData,M> > >& celldata, double)
    {
        const int n = _pos.size();
        celldata.resize(n);
        for(int i=0;i<n;++i) {
            celldata[i].reset(new CellData<NData,M>(_pos[i]));
        }
        dbg<<"ngal = "<<celldata.size()<<std::endl;
    }

    template <int M>
    void buildCellData(double minsep, double binsize,
                       std::vector<boost::shared_ptr<CellData<EData,M> > >& celldata, double& vare)
    {
        Assert(_shear.size() == _pos.size());
        Assert(_weight.size() == 0 || _weight.size() == _pos.size());
        const int n = _pos.size();
        celldata.resize(n);
        double sumw=0.;
        vare=0.;
        for(int i=0;i<n;++i) {
            double w = _weight.size() > 0 ? _weight[i] : 1.;
            celldata[i].reset(new CellData<EData,M>(_pos[i],_shear[i],w));
            sumw += w;
            vare += w*w*norm(_shear[i]);
        }
        vare /= 2*sumw*sumw/celldata.size(); // 2 because want var per component
        dbg<<"ngal = "<<celldata.size()<<", totw = "<<sumw<<std::endl;
    }

    bool useRaDec() const { return _use_ra_dec; }

private:

    bool _use_ra_dec;
    std::vector<Position<Flat> > _pos;
    std::vector<std::complex<double> > _shear;
    std::vector<double> _weight;
};

void ReadInputFiles(std::vector<boost::shared_ptr<InputFile> >& files, const ConfigFile& params)
{
    std::vector<std::string> filename = params["file_name"];
    const int nfiles = filename.size();
    for (int i=0; i<nfiles; ++i) {
        boost::shared_ptr<InputFile> file(new InputFile(filename[i],params));
        files.push_back(file);
    }
}

#endif
